rm(list = ls())

#setwd("~/Dropbox/Energiewende_Project/All_Replication_Files")
setwd("C:/Users/chris/Dropbox/Energiewende_Project/All_Replication_Files")

#Loading in Packages ----------------------------------------------------------------------------------

library(dplyr)
library(sandwich)
library(lmtest)
library(margins)
library(ggplot2)
library(scales)
library(gridExtra)
library(readr)
library(stargazer)
library(reshape2) 
library(stringr)
library(cregg)
library(kableExtra)
library(knitr)
library(modelsummary)
options(knitr.table.format = "latex")

####################################################################################################################################
############################################Coal Exit Experiment####################################################################
####################################################################################################################################
###########Prep Data:####
data <- read_csv("data.coal.cleaned.final.csv", col_types = cols(.default = "?"))

data$coalexit <- relevel(as.factor(data$coalexit), "Control")



###########Helper Functions:#####
#main regression:
reg_coal <- function(data, subgroup){
  
  datac <- data[data[,subgroup] == 1,]
  fmla= function(y, x=x.coal) as.formula(paste0(y, '~', x))
  mods.coal = lapply(y.variables, fmla)
  
  mod.fit = function(fmla) lm(fmla, data=datac)
  
  fits.coal = lapply(mods.coal, mod.fit)
  
  return(fits.coal)
}

#with covariates

reg_coal_cov <- function(data, subgroup){
  
  datac <- data[data[,subgroup] == 1,]
  fmla= function(y, x=x.coal) as.formula(paste0(y, '~', x, '+ region + age + gender + incomequintiles'))
  mods.coal = lapply(y.variables, fmla)
  
  mod.fit = function(fmla) lm(fmla, data=datac)
  
  fits.coal = lapply(mods.coal, mod.fit)
  
  return(fits.coal)
}

#with heterogeneity interaction terms:


reg_coal_int<- function(data,inter, subgroup){
  
  datac <- data[data[,subgroup] == 1,]
  fmla= function(y, x=x.coal) as.formula(paste0(y, '~', x, '+',inter, '+',inter,'*',x))
  mods.coal = lapply(y.variables, fmla)
  
  mod.fit = function(fmla) lm(fmla, data=datac)
  
  fits.coal = lapply(mods.coal, mod.fit)
  
  return(fits.coal)
}

reg_coal_int_cov<- function(data,inter, subgroup){
  
  datac <- data[data[,subgroup] == 1,]
  fmla= function(y, x=x.coal) as.formula(paste0(y, '~', x, '+',inter, '+',inter,'*',x,'+ region + age + gender + incomequintiles'))
  mods.coal = lapply(y.variables, fmla)
  
  mod.fit = function(fmla) lm(fmla, data=datac)
  
  fits.coal = lapply(mods.coal, mod.fit)
  
  return(fits.coal)
}


##robustness to other attitudes


reg_coal_int_rob <- function(data,inter, subgroup){
  
  datac <- data[data[,subgroup] == 1,]
  fmla= function(y, x=x.coal) as.formula(paste0(y, '~', x, '+',inter, '+',inter,'*',x, "+ climatehumans*",x,"+climate_policyn*",x,"+ left_right_4*",x,"+ party_2021_2*",x))
  mods.coal = lapply(y.variables, fmla)
  
  mod.fit = function(fmla) lm(fmla, data=datac)
  
  fits.coal = lapply(mods.coal, mod.fit)
  
  return(fits.coal)
}

#data prep for plotting:
data.prep <- function(reg){
  regsum <- summary(reg)
  coef <- as.matrix(regsum$coefficients[2:5,1:2])  
  coef <- as.data.frame(cbind(c("Compensation for Consumers Treatment", "Compensation for Investors Treatment", 
                                "Compensation for Regions Treatment", "Compensation for Workers Treatment"), coef))
  coef$Estimate <- as.numeric(paste(coef$Estimate))
  coef$`Std. Error`<- as.numeric(paste(coef$`Std. Error`))
  coef$confupper <- coef$Estimate + coef$`Std. Error`*1.96
  coef$conflower <- coef$Estimate - coef$`Std. Error`*1.96
  coef$V1 <- as.factor(coef$V1)
  coef$V1 <- factor(coef$V1,
                    levels = levels(coef$V1)[c(1,2,3,4,5)])
  return(coef)
}

#regression tables:

#plot results:
plot_basic <- function(regout, fileadd, dvadd){
  data.temp <- data.prep(regout)
  
  pl <-  ggplot(data.temp)+
    geom_point(aes(x = Estimate, y = V1))+
    geom_segment(aes(x = conflower, xend = confupper,
                     y = V1, yend = V1))+
    geom_vline(xintercept = 0, lty = "dotted")+
    theme_bw()+
    ggtitle(paste0("Treatment Effects Coal Vignette, \n", dvadd))
  
  
  pl <- pl + theme(text = element_text(size=16),
                   plot.title = element_text(size=18, hjust = 0.5),
                   axis.text=element_text(size=16),
                   axis.title = element_text(size=16))+
    scale_y_discrete(labels = function(x) str_wrap(x, width = 15))+
    xlab("Average Treatment Effects")+
    ylab("")
  
  print(pl)
  
  ggsave(filename = paste0("Graphs/coal_vignette.", fileadd, ".png"),
         plot = pl,
         width = 12, height = 8)
}


#heterogeneity analysis:
drawplots_coal <- function(res, n.res, colors, legendtitle ,addtitle=NULL, fileadd = NULL){
  for(k in 1:4){
    plot.datac <- lapply(res,"[[",k)
    
    plot.datap <- lapply(plot.datac, data.prep)
    
    data.temp <- as.data.frame(plot.datap[[1]])
    
    pl <-  ggplot(data.temp)+
      geom_point(aes(x = Estimate, y = V1, color = paste0(colors[1] , ", n= ", n.res[[1]]), shape = paste0(colors[1] , ", n= ", n.res[[1]])))+
      geom_segment(aes(x = conflower, xend = confupper,
                       y = V1, yend = V1, color = paste0(colors[1] , ", n= ", n.res[[1]])))+
      geom_vline(xintercept = 0, lty = "dotted")+
      theme_bw()+
      ggtitle(paste0("Treatment Effects Coal Vignette", addtitle, "\n ", dependent.variables[k]))
    
    
    
    
    if(length(res) > 1){
      data.temp <- as.data.frame(plot.datap[[2]])
      pl <- pl +  geom_point(data = data.temp, aes(x = Estimate, y = V1, color = paste0(colors[2] , ", n= ", n.res[[2]]), shape = paste0(colors[2] , ", n= ", n.res[[2]])), position = position_nudge(y = -0.1))+
        geom_segment(data = data.temp,aes(x = conflower, xend = confupper,
                                          y = V1, yend = V1, color = paste0(colors[2] , ", n= ", n.res[[2]])), position = position_nudge(y = -0.1))
    }
    
    if(length(res) > 2){
      data.temp <- as.data.frame(plot.datap[[3]])
      pl <- pl +  geom_point(data = data.temp, aes(x = Estimate, y = V1, color = paste0(colors[3] , ", n= ", n.res[[3]])), position = position_nudge(y = +0.1))+
        geom_segment(data = data.temp,aes(x = conflower, xend = confupper,
                                          y = V1, yend = V1, color = paste0(colors[3] , ", n= ", n.res[[3]])), position = position_nudge(y = +0.1))
    }
    
    if(length(res) > 3){
      data.temp <- as.data.frame(plot.datap[[4]])
      pl <- pl +  geom_point(data = data.temp, aes(x = Estimate, y = V1, color = paste0(colors[4] , ", n= ", n.res[[4]])), position = position_nudge(y = -0.2))+
        geom_segment(data = data.temp,aes(x = conflower, xend = confupper,
                                          y = V1, yend = V1, color = paste0(colors[4] , ", n= ", n.res[[4]])), position = position_nudge(y = -0.2))
    }
    
    
    if(length(res) > 4){
      data.temp <- as.data.frame(plot.datap[[5]])
      pl <- pl +  geom_point(data = data.temp, aes(x = Estimate, y = V1, color = paste0(colors[5] , ", n= ", n.res[[5]])), position = position_nudge(y = +0.2))+
        geom_segment(data = data.temp,aes(x = conflower, xend = confupper,
                                          y = V1, yend = V1, color = paste0(colors[5] , ", n= ", n.res[[5]])), position = position_nudge(y = +0.2))
    }
    
    if(length(res) > 5){
      data.temp <- as.data.frame(plot.datap[[6]])
      pl <- pl +  geom_point(data = data.temp, aes(x = Estimate, y = V1, color = paste0(colors[6] , ", n= ", n.res[[6]])), position = position_nudge(y = -0.3))+
        geom_segment(data = data.temp,aes(x = conflower, xend = confupper,
                                          y = V1, yend = V1, color = paste0(colors[6] , ", n= ", n.res[[6]])), position = position_nudge(y = -0.3))
    }
    
    pl <- pl + labs(color='Group') + theme(text = element_text(size=16),
                                           plot.title = element_text(size=18, hjust = 0.5),
                                           axis.text=element_text(size=16),
                                           axis.title = element_text(size=16),
                                           legend.position = "bottom",
                                           legend.text = element_text(size=16))+
      scale_y_discrete(labels = function(x) str_wrap(x, width = 15))+
      xlab("Average Treatment Effects")+
      ylab("") + 
      guides(color=guide_legend(title=legendtitle), 
             shape = guide_legend(title=legendtitle))+
      scale_color_grey(start = 0.1, end = 0.5)
    
    print(pl)
    
    ggsave(filename = paste("Graphs/coal_vignette.", fileadd, file.names.y[k],".png", sep=""),
           plot = pl,
           width = 12, height = 8)
    
  }
  
}


###########Define Terms and Vectors:#####
y.variables <- c("coalexitpremax", "coalexitsooner", "E3_exitsupport","coalexitsupport")

x.coal = 'coalexit'

dependent.variables <- c("2050 - Preferred Year of Exit", "Coal Exit Sooner than 2038", "Support for Full Exit, Continuous", "Support for Full Exit, Binary")

file.names.y <- c("maxyear", "sooner", "supportc", "supportb")

supportinteractions <- c("highsupportinterventionecon", "highsupportinterventionfrail", "highsupportinterventionworkers")


###########Basic Regression:####


fits.coal2 <- reg_coal_cov(data, "all")
#drawing plots of four main regressions:

for(i in 1:4){
  plot_basic(fits.coal2[[i]], file.names.y[i], dependent.variables[i])
}




####heterogeneity plots with split samples

plot.data <- list()
n.subsets <- list()
subsets <-  c("highsupportintervention", "lowsupportintervention")

for(i in 1:length(subsets)){
  regressionsout <- reg_coal_cov(data, subsets[i])
  plot.data[[i]] <- regressionsout
  n.subsets[[i]] <- length(regressionsout[[1]]$fitted.values)
}

drawplots_coal(plot.data, n.subsets, c("High", "Low") , paste0("Belief State Intervention"),  " by Belief in State Intervention", "govtint") 


##########Creating Regression Tables for Paper:###############

# #table 1: baseline and interaction with belief in state management

for(i in 1:4){
  reg.1 <- reg_coal(data, "all")[[i]]

  reg.2 <-  reg_coal_cov(data, "all")[[i]]

  reg.3 <- reg_coal_int(data, "highsupportintervention", "all")[[i]]

  reg.4 <- reg_coal_int_cov(data, "highsupportintervention", "all")[[i]]

  reg.5 <- reg_coal_int(data, "supportintervention", "all")[[i]]

  reg.6 <- reg_coal_int_cov(data, "supportintervention", "all")[[i]]


  stargazer(reg.2, reg.4, reg.6,
            keep=c("coalexitCompensation for Consumers",
                   "coalexitCompensation for Investors",
                   "coalexitCompensation for Regions",
                   "coalexitCompensation for Workers",
                   "highsupportintervention",
                   "coalexitCompensation for Consumers:highsupportintervention",
                   "coalexitCompensation for Investors:highsupportintervention",
                   "coalexitCompensation for Regions:highsupportintervention",
                   "coalexitCompensation for Workers:highsupportintervention",
                   "supportintervention",
                   "coalexitCompensation for Consumers:supportintervention",
                   "coalexitCompensation for Investors:supportintervention",
                   "coalexitCompensation for Regions:supportintervention",
                   "coalexitCompensation for Workers:supportintervention"),
            covariate.labels =c("Compensation for Consumers","Compensation for Investors",
                                 "Compensation for Regions","Compensation for Workers",
                                 "Belief State Int. Bin.", "Belief State Int. Cont.",
                                "Comp. Consumers * Belief State Int. Bin.",
                                 "Comp. Investors * Belief State Int. Bin.",
                                 "Comp. Regions * Belief State Int. Bin.",
                                 "Comp. Workers * Belief State Int. Bin.",
                                  "Comp. Consumers * Belief State Int. Cont.",
                                 "Comp. Investors * Belief State Int. Cont.",
                                 "Comp. Regions * Belief State Int. Cont.",
                                 "Comp. Workers * Belief State Int. Cont."),
            add.lines = list(c("Demographic Controls",  "Yes",  "Yes",  "Yes")),
            font.size="tiny",
            dep.var.labels = dependent.variables[i],
            out=paste0("Tables/coal_vignette,",file.names.y[i],".tex")
  )


}

##table 1 repeated for only those who passed the mock vignette

for(i in 1:4){
  reg.1.mv <- reg_coal(data, "coalattention")[[i]]

  reg.2.mv <-  reg_coal_cov(data, "coalattention")[[i]]

  reg.3.mv <- reg_coal_int(data, "highsupportintervention", "coalattention")[[i]]

  reg.4.mv <- reg_coal_int_cov(data, "highsupportintervention", "coalattention")[[i]]

  reg.5.mv <- reg_coal_int(data, "supportintervention", "coalattention")[[i]]

  reg.6.mv <- reg_coal_int_cov(data, "supportintervention", "coalattention")[[i]]


  stargazer(reg.2.mv,  reg.4.mv,  reg.6.mv,
            keep=c("coalexitCompensation for Consumers",
                   "coalexitCompensation for Investors",
                   "coalexitCompensation for Regions",
                   "coalexitCompensation for Workers",
                   "highsupportintervention",
                   "coalexitCompensation for Consumers:highsupportintervention",
                   "coalexitCompensation for Investors:highsupportintervention",
                   "coalexitCompensation for Regions:highsupportintervention",
                   "coalexitCompensation for Workers:highsupportintervention",
                   "supportintervention",
                   "coalexitCompensation for Consumers:supportintervention",
                   "coalexitCompensation for Investors:supportintervention",
                   "coalexitCompensation for Regions:supportintervention",
                   "coalexitCompensation for Workers:supportintervention"),
            covariate.labels =c("Compensation for Consumers","Compensation for Investors",
                                "Compensation for Regions","Compensation for Workers",
                                "Belief State Int. Bin.",  "Belief State Int. Cont.","Comp. Consumers * Belief State Int. Bin.",
                                "Comp. Investors * Belief State Int. Bin.",
                                "Comp. Regions * Belief State Int. Bin.",
                                "Comp. Workers * Belief State Int. Bin.",
                                "Comp. Consumers * Belief State Int. Cont.",
                                "Comp. Investors * Belief State Int. Cont.",
                                "Comp. Regions * Belief State Int. Cont.",
                                "Comp. Workers * Belief State Int. Cont."),
            add.lines = list(c("Demographic Controls", "Yes", "Yes", "Yes")),
            font.size="tiny",
            dep.var.labels = dependent.variables[i],
            out=paste0("Tables/coal_vignette_mv_passed,",file.names.y[i],".tex")
  )


}

#split sample tables

list.tablesplit.h <- reg_coal_cov(data, "highsupportintervention")
list.tablesplit.l <- reg_coal_cov(data, "lowsupportintervention")

list.tablesplit.h.a <- reg_coal_cov(data[data$coalattention==1,], "highsupportintervention")
list.tablesplit.l.a <- reg_coal_cov(data[data$coalattention==1,], "lowsupportintervention")

for(i in 1:4){
  
  table_split <- modelsummary(list(list.tablesplit.h[[i]], list.tablesplit.l[[i]], list.tablesplit.h.a[[i]], list.tablesplit.l.a[[i]]),
                              stars=T,
                              coef_rename = c("coalexitCompensation for Consumers Treatment" = "Compensation for Consumers",
                                              "coalexitCompensation for Investors Treatment" = "Compensation for Investors",
                                              "coalexitCompensation for Regions" = "Compensation for Regions",
                                              "coalexitCompensation for Workers" = "Compensation for Workers"),
                              coef_omit = c("region|age|gender|income"),
                              output="latex")
  
  header.t1 <- paste("DV:", dependent.variables[i])
  
  tableheader.ts1 <- c(c(" " = 1,  header.t1 = 4 ))
  
  # set vector names 
  names(tableheader.ts1) <- c(" ", header.t1)
  
  high_a <- "Passed Attention"
  all <- "All"
  
  tableheader.ts2 <- c(" " = 1, all=2, high_a = 2)
  names(tableheader.ts2) <- c(" ", all, high_a)
  
  high_b <- "High Belief St. Int."
  low_b <- "Low Belief St. Int."
  
  tableheader.ts3 <- c(" " = 1, high_b =1, low_b = 1, high_b=1, low_b=1)
  names(tableheader.ts3) <- c(" ", high_b, low_b, high_b, low_b)
  
  regtable.ts1 <- table_split %>% kableExtra::add_header_above(header=tableheader.ts3)%>% kableExtra::add_header_above(header=tableheader.ts2) %>% kableExtra::add_header_above(header=tableheader.ts1) 
  
  
  save_kable(regtable.ts1, paste0("Tables/coal_vignette_split_",file.names.y[i],".tex"))
  
}



###Table 3: Robust interactions with interactions with other attitude variables


for(i in 1:4){
  
  
  header.t1 <- paste("DV:", dependent.variables[i])
  
  tableheader.tr1 <- c(c(" " = 1,  header.t1 = 1))
  
  reg.1.r <- reg_coal_int_rob(data, "highsupportintervention", "all")[[i]]
  
  
  
  table.rf1 <- modelsummary(list(reg.1.r), stars=T,
                            gof_omit="AIC|BIC",
                            coef_rename = c("coalexitCompensation for Consumers" = "Compensation for Consumers",
                                            "coalexitCompensation for Investors" = "Compensation for Investors",
                                            "coalexitCompensation for Regions" = "Compensation for Regions",
                                            "coalexitCompensation for Workers" = "Compensation for Workers",
                                            "party_2021_2" = "Party Vote 2021",
                                            "left_right_4" = "Left-Right Leaning",
                                            "climate_policyn"= "Climate Policy Support",
                                            "climatehumans" = "Belief in MM CC",
                                            "highintervention" = "Belief in St. Int. Bin.",
                                            "supportintervention" = "Belief in St. Int. Cont.",
                                            "region" = "Region "),
                            estimate  = "{estimate} ({std.error})",
                            statistic = NULL,
                            output="latex")
  
  
  
  regtable.rf1 <- table.rf1 %>% kableExtra::add_header_above(header=tableheader.tr1) 
  
  
  save_kable(regtable.rf1, paste0("Tables/coal_vignette_robust1,",file.names.y[i],".tex"))
  
  
  
  reg.2.r <- reg_coal_int_rob(data, "supportintervention", "all")[[i]]
  
  table.rf2 <- modelsummary(list(reg.2.r), stars=T,
                            gof_omit="AIC|BIC",
                            coef_rename = c("coalexitCompensation for Consumers" = "Compensation for Consumers",
                                            "coalexitCompensation for Investors" = "Compensation for Investors",
                                            "coalexitCompensation for Regions" = "Compensation for Regions",
                                            "coalexitCompensation for Workers" = "Compensation for Workers",
                                            "party_2021_2" = "Party Vote 2021",
                                            "left_right_4" = "Left-Right Leaning",
                                            "climate_policyn"= "Climate Policy Support",
                                            "climatehumans" = "Belief in MM CC",
                                            "highintervention" = "Belief in St. Int. Bin.",
                                            "supportintervention" = "Belief in St. Int. Cont.",
                                            "region" = "Region "),
                            estimate  = "{estimate} ({std.error})",
                            statistic = NULL,
                            output="latex")
  
  
  
  regtable.rf2 <- table.rf2 %>% kableExtra::add_header_above(header=tableheader.tr1) 
  
  
  save_kable(regtable.rf2, paste0("Tables/coal_vignette_robust2,",file.names.y[i],".tex"))
  

  
}

#checking just industry intervention


reg.1.i <- reg_coal(data,  "ind_help_supp")

reg.2.i <- reg_coal_cov(data,  "ind_help_supp")

modelsummary(list(reg.1.i[[4]],reg.2.i[[4]]), stars=T, output="latex")


